bathymap-for-amalga

24 Feb 2022

# library(tidyverse)
# library(ncdf4)
# library(raster)

# Let's see if this DEM from NOAA is useful
# juneau <- raster("JuneauBathy.tif")

# plot(juneau)

Filter out the land values by setting >=0

leaning on this: https://www.benjaminbell.co.uk/2019/08/bathymetric-maps-in-r-colour-palettes.html

To figure out a color scheme, need min and max

# Calculate min and max values
# mi <- cellStats(juneau, stat="min")-100
# ma <- cellStats(juneau, stat="max")+100

    # # Break points sequence for below sea level
    # s1 <- seq(from=mi, to=0, by=0 - mi / 50)
    # # Break points sequence for above sea level
    # s2 <- seq(from=0, to=ma, by=ma / 50)
    # 
    #     # Round sequence to nearest 100
    # s1 <- round(s1, -1)
    # s2 <- round(s2, -1)
    # 
    #     # Only show unique numbers
    # s1 <- unique(s1)
    # s2 <- unique(s2)
    # 
    # s1
    # 
    # s2

    # Combine sequences and remove the first value from second sequence
    # s3 <- c(s1, s2[-1])

    # Plot
    # plot(juneau, col=c(blue.col(50), terrain.colors(50)), breaks=s3)

# crop <- extent(-134.8,-134.6,58.45,58.5)
# 
# extent(juneau)
# 
# june.zoom <- setExtent(juneau, crop)
# 
# plot(june.zoom, col=c(blue.col(50), terrain.colors(50)), breaks=s3)

in ggplot

https://stackoverflow.com/questions/33227182/how-to-set-use-ggplot2-to-map-a-raster

library(rnaturalearth)

Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`

Attaching package: ‘rnaturalearth’

The following object is masked from ‘package:rnaturalearthdata’:

    countries110

datafold <- "../data/JuneauBathy.tif"
test <- raster(datafold) 

test_spdf <- as(test, "SpatialPixelsDataFrame")
test_df <- as.data.frame(test_spdf)
colnames(test_df) <- c("value", "x", "y")

zoom.df <- test_df %>%
  filter(x > -134.95, x < -134.76, y > 58.4, y < 58.55)
  #filter(x > -135, x < -134, y > 58.2, y < 58.7) 

zoom.df.w.gray <- zoom.df %>%
  mutate(value = ifelse(value >=0, NA, value)) # set land to NA to shade it gray

# fully zoomed out
full_size_raster <- test_df %>%
  mutate(value = ifelse(value >=0, NA, value)) # set land to NA to shade it gray

# minimum value below sea level
mi <- min(zoom.df$value)

# Break points sequence for below sea level
s1 <- seq(from=mi, to=0, by=0 - mi / 50)

depth.scale <- round(s1, 0)
depth.scale <- unique(depth.scale)

depth.scale

 [1] -262 -257 -252 -247 -241 -236 -231 -226 -220 -215 -210 -205 -199 -194 -189 -184 -178 -173 -168 -163 -157 -152 -147 -142 -136 -131 -126 -121
[29] -115 -110 -105 -100  -94  -89  -84  -79  -73  -68  -63  -58  -52  -47  -42  -37  -31  -26  -21  -16  -10   -5    0

breaks = levels(depth.scale)[floor(seq(1, nlevels(depth.scale), length.out = 10))]

new.depth.scale <- depth.scale[seq(1, length(depth.scale), 3)]

zoomed.out <- ggplot() +  
  geom_tile(data=zoom.df.w.gray, aes(x=x, y=y, fill=value), alpha=0.8) + 
  scale_fill_steps2(low = "navyblue",
                    mid = "dodgerblue4",
  high = "lightsteelblue1",
  midpoint = -155,
  na.value = "grey50",
  breaks = new.depth.scale) +
  coord_equal() +
  theme(panel.border = element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 12),
        axis.title.x = element_text(size = 14, margin = margin(t=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10)),
        axis.line = element_line(color = "black"),
        legend.position="right") +
      theme(legend.text = element_text(size = 8)) + 
      theme(legend.key.height=unit(1.8, "cm")) +
      labs(fill = "Depth (m)",
         x = "Longitude (W)",
         y = "Latitude (N)")

zoomed.out

#ggsave("pdf_outputs/amalgaBathy.pdf")

That’s more zoomed out, but gives a better perspective…

Then here’s more zoomed in on Amalga:

more.zoom.df <- test_df %>%
  filter(x > -134.85, x < -134.77, y > 58.43, y < 58.51) %>%
  mutate(value = ifelse(value >=0, NA, value)) # set land to NA to shade it gray

zoomed.map <- ggplot() +  
  geom_tile(data=more.zoom.df, aes(x=x, y=y, fill=value), alpha=0.8) + 
  scale_fill_steps2(low = "navyblue",
                    mid = "dodgerblue4",
  high = "lightsteelblue1",
  midpoint = -97,
  #space = "Lab",
  na.value = "grey50",
  breaks = new.depth.scale) +
  coord_equal() +
  coord_equal() +
  theme(panel.border = element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 12),
        axis.title.x = element_text(size = 14, margin = margin(t=10)),
        axis.title.y = element_text(size = 14, margin = margin(r=10)),
        axis.line = element_line(color = "black"),
        legend.position="right") +
      theme(legend.text = element_text(size = 8)) + 
      theme(legend.key.height=unit(1.8, "cm")) +
      labs(fill = "Depth (m)",
         x = "Longitude (W)",
         y = "Latitude (N)")

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

zoomed.map

#ggsave("amalgaBathyZoom.pdf")

Add the transect line to this map?

# read in the data
# this dataframe is produced in the CTD analysis on the VM
# 08-ctd-cast-data.Rmd
ctd <- read_csv("../data/ctdDataframe.csv")

Rows: 2538 Columns: 16── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr   (2): ctd_sample, tide
dbl  (13): lat, long, depth_m, salinity, density, pressure_decibar, temp_c, conductivity, sp_conduct, sound_velocity, duration, id, distance
dttm  (1): time
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# just one of the transects (not both AM and PM)
am <- ctd %>%
  filter(tide == "AM_incoming") %>%
  dplyr::select(lat, long, id, distance) %>%
  unique() %>%
  filter(id != 1)

am %>%
  ggplot(aes(x = long, y = lat, label = id)) +
  geom_text()

zoomed.map +
  geom_point(data = am, aes(x = long, y = lat, label = id), size = 1, color = "firebrick2")

Warning: Ignoring unknown aesthetics: label

#ggsave("pdf_outputs/amalgaTransectBathy.pdf", height = 6, width = 7)

Zoom more

highest.res <- more.zoom.df %>%
  filter(x > -134.83, x < -134.779, y > 58.47, y < 58.51) %>%
  mutate(value = ifelse(value >=0, NA, value)) # set land to NA to shade it gray

highest.zoom.plot <- ggplot() +  
  geom_tile(data=highest.res, aes(x=x, y=y, fill = value, color = value)) + 
  scale_fill_steps2(low = "navyblue",
                    mid = "dodgerblue4",
  high = "lightsteelblue1",
  midpoint = -99,
  space = "Lab",
  na.value = "grey50",
  breaks = new.depth.scale) +
  scale_color_steps2(low = "navyblue",
                    mid = "dodgerblue4",
  high = "lightsteelblue1",
  midpoint = -99,
  na.value = "grey50",
  breaks = new.depth.scale) +
  coord_equal() +
  theme(panel.border = element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 8),
        axis.title.x = element_text(size = 10, margin = margin(t=10)),
        axis.title.y = element_text(size = 10, margin = margin(r=10)),
        axis.line = element_line(color = "gray50"),
        legend.position="right") +
       theme(legend.text = element_text(size = 6)) + 
      # theme(legend.key.height=unit(1, "cm")) +
      labs(fill = "Depth (m)",
         x = "Longitude (W)",
         y = "Latitude (N)") +
  guides(color = F) +
  scale_y_continuous(expand = c(0,0), breaks = c(58.48, 58.49, 58.50)) +
  scale_x_continuous(expand = c(0,0))
  

highest.zoom.plot

transect2021 <- highest.zoom.plot +
  geom_point(data = am, aes(x = long, y = lat), size = 1, color = "red") +
  #annotate(geom = "label", x = -134.795, y = 58.504, label = "Amalga Harbor", size = 4) +
  #annotate(geom = "label", x = -134.81, y = 58.494, label = "2021 transect", size = 3, fontface = "bold") +
  annotate(geom = "label", x = -134.78999, y = 58.496, label = "pens", size = 3, color = "black") +
    annotate(geom = "label", x = -134.808, y = 58.492, label = "1000m", size = 3, color = "black") +
      annotate(geom = "label", x = -134.825, y = 58.485, label = "2000m", size = 3, color = "black") +
  scale_x_continuous(limits = c(-134.83, -134.78), breaks = c(-134.82, -134.80, -134.78)) +
    scale_y_continuous(limits = c(58.47, 58.51)) 

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

  #labs(title = "2021 transect") +
  #scale_fill_continuous(breaks = c(0, -10, -25, -50, -75, -100, -125, -150, -200))
  
  
#ggsave("pdf_outputs/amalgaTransectBathyZoom.pdf", width = 6, height = 5)

at 560 m, it’s basically the edge of Amalga Harbor.

Add the transect to the zoomed out map:

map2021_out <- zoomed.out +
  geom_point(data = am, aes(x = long, y = lat), size = 1, color = "black") +
  annotate(geom = "text", x = -134.786, y = 58.51, label = "Amalga Harbor", size = 4, color = "white", fontface = "bold") +
  annotate(geom = "label", x = -134.83, y = 58.494, label = "2021 transect", size = 3, fontface = "bold")
  

With these maps, I should be able to combine the two transect years with cowplot or similar.

For the 2022 data, I need the equivalent from the CTD casts for the transect.

CTD data from 2022

ctd.2022 <- read_csv("../data/2022ctdDataframe.csv")

Rows: 1360 Columns: 15── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr   (2): ctd_sample, tide
dbl  (12): lat, long, depth_m, salinity, density, pressure_decibar, temp_c, conductivity, sp_conduct, sound_velocity, duration, depth
dttm  (1): time
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ctd.2022

# grab the data for May 5 and one tide for data points for the map
one.set.2022 <- ctd.2022 %>%
  filter(str_detect(time, "2022-05-05") &
    tide == "AM_outgoing") %>%
  dplyr::select(ctd_sample, lat, long) %>%
  unique() %>%
  filter(!ctd_sample %in% c("163024", "191935", "175337", "185047")) %>%
  mutate(Year = "2022")

fp <- am %>%
  mutate(Year = "2021") %>%
  bind_rows(one.set.2022)


amalga_plot <- highest.zoom.plot +
  new_scale_color() +
  geom_point(data = fp, aes(x = long, y = lat, shape = Year, color = Year), size = 1.5) +
  annotate(geom = "label", x = -134.789, y = 58.496, label = "pens", size = 3, color = "black") +
    annotate(geom = "label", x = -134.808, y = 58.4865, label = "1000m", size = 3, color = "black") +
      annotate(geom = "label", x = -134.826, y = 58.485, label = "2000m", size = 3, color = "black") +
  scale_color_manual(values = c("coral3", "gold")) +
  scale_shape_manual(values = c(16,17)) +
  theme(
    legend.position = "bottom",
    legend.spacing.x = unit(0, "cm"),
    legend.key = element_rect(fill = "white"),
    legend.spacing = unit(3, "cm")
  ) +
  guides(fill = guide_legend(label.position = "bottom",
                             nrow = 1,
                             keywidth = unit(0.5, "cm"),
                             title.hjust = unit(0.2, "cm"),
                             title.vjust = unit(0.4, "cm")),
         color = guide_legend(override.aes = list(size = 5)))
                             

amalga_plot

ggsave("pdf_outputs/amalgaTransectBathyZoom.png", width = 6, height = 5)

# z.o <- zoomed.out +
#   geom_point(data = one.set.2022, aes(x = long, y = lat), size = 1, color = "white") +
#   annotate(geom = "label", x = -134.8, y = 58.51, label = "Amalga Harbor", size = 4, color = "black") +
#   annotate(geom = "label", x = -134.86, y = 58.48, label = "background", size = 3, color = "black") +
#     annotate(geom = "label", x = -134.80, y = 58.47, label = "transects", size = 3, color = "black") +
#   theme(
#      legend.position = "none"
#   )
# 
#   
# z.i.22 <- highest.zoom.plot +
#   geom_point(data = one.set.2022, aes(x = long, y = lat), size = 1, color = "red") +
#   annotate(geom = "label", x = -134.78999, y = 58.496, label = "pens", size = 3, color = "black") +
#     annotate(geom = "label", x = -134.808, y = 58.492, label = "1000m", size = 3, color = "black") +
#       annotate(geom = "label", x = -134.825, y = 58.485, label = "2000m", size = 3, color = "black") +
#   scale_x_continuous(limits = c(-134.83, -134.78), breaks = c(-134.82, -134.80, -134.78)) +
#     scale_y_continuous(limits = c(58.47, 58.51)) +
#   theme(
#   #   legend.position = "none",
#   #   axis.text.x = element_blank(),
#   #   axis.title.x = element_blank(),
#   #   axis.ticks.x = element_blank(),
#     axis.title.y = element_blank()  
#     )
# 
# z.i.21 <- highest.zoom.plot +
#   geom_point(data = am, aes(x = long, y = lat), size = 1, color = "red") +
#   annotate(geom = "label", x = -134.78999, y = 58.496, label = "pens", size = 3, color = "black") +
#     annotate(geom = "label", x = -134.808, y = 58.492, label = "1000m", size = 3, color = "black") +
#       annotate(geom = "label", x = -134.825, y = 58.485, label = "2000m", size = 3, color = "black") +
#   scale_x_continuous(limits = c(-134.83, -134.78), breaks = c(-134.82, -134.80, -134.78)) +
#     scale_y_continuous(limits = c(58.47, 58.51)) +
#    theme(
#     #legend.position = "none",
#     axis.text.x = element_blank(),
#     axis.title.x = element_blank(),
#     axis.ticks.x = element_blank(),
#     axis.title.y = element_blank()
#   )

library(ggOceanMaps)

ggOceanMaps: Setting data download folder to a temporary folder /var/folders/_5/zm9fzz0d1j139pfm56sqzlfs8qyj5v/T//RtmpQYxlZb. This
means that any downloaded map data need to be downloaded again when you restart R. To avoid this problem, change the default path
to a permanent folder on your computer. Add following lines to your .Rprofile file: {.ggOceanMapsenv <- new.env();
.ggOceanMapsenv$datapath <- 'YourCustomPath'}. You can use usethis::edit_r_profile() to edit the file. '~/ggOceanMapsLargeData'
would make it in a writable folder on most operating systems.

Attaching package: ‘ggOceanMaps’

The following object is masked from ‘package:ggthemes’:

    theme_map

Combine the ak plot and the Amalga plot

(ak_plot | interm_p) / amalga_plot2 + plot_layout(nrow = 2, heights = c(2,3)) +
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(size = 8),
        plot.margin = margin(0.1, 0.1, 0.1, 0.1, "cm"))

ggsave("pdf_outputs/amalga_combined_maps.png", width = 7, height = 7)

NOT IN USE BELOW THIS LINE

Bathy transect

library(marmap)

# load datasets
#   data(nw.atlantic); as.bathy(nw.atlantic) -> atl
#   data(nw.atlantic.coast)
# 
# # Example 1. get.transect(), without use of locator()
#   get.transect(atl, -65, 43,-59,40) -> test ; plot(test[,3]~test[,2],type="l")
#   get.transect(atl, -65, 43,-59,40, distance=TRUE) -> test ; plot(test[,4]~test[,3],type="l")
# 
# # Example 2. get.transect(), without use of locator(); pretty plot
#   par(mfrow=c(2,1),mai=c(1.2, 1, 0.1, 0.1))
#   plot(atl, deep=-6000, shallow=-10, step=1000, lwd=0.5, col="grey50",drawlabels=TRUE)
#   lines(nw.atlantic.coast)
# 
#   get.transect(atl, -75, 44,-46,32, loc=FALSE, dis=TRUE) -> test
#   points(test$lon,test$lat,type="l",col="blue",lwd=2,lty=2)
#   plotProfile(test)
# 
# # Example 3. get.transect(), with use of locator(); pretty plot
# ## Not run: 
#   par(mfrow=c(2,1),mai=c(1.2, 1, 0.1, 0.1))
#   plot(atl, deep=-6000, shallow=-10, step=1000, lwd=0.5, col="grey50",drawlabels=TRUE)
#   lines(nw.atlantic.coast)
#   
#   get.transect(atl, loc=TRUE, dis=TRUE, col=2, lty=2) -> test
#   plotProfile(test)
    
## End(Not run)

Very cool, so if I can get my data into that format, I should be able to do something similar.

# -134.7925, -134.8213, 58.4946, 58.4848
# load dataset
juneau <- getNOAA.bathy(lon1 = -140, lon2 = -130,
lat1 = 60, lat2 = 55, resolution = 1)

plot(juneau)
summary(juneau)

# Creating a custom palette of blues
blues <- c("lightsteelblue4", "lightsteelblue3",
"lightsteelblue2", "lightsteelblue1")
# Plotting the bathymetry with different colors for land and sea
plot(juneau, image = TRUE, land = TRUE, lwd = 0.1,
bpal = list(c(0, max(juneau), "grey"),
c(min(juneau),0,blues)))
# Making the coastline more visible
plot(juneau, deep = 0, shallow = 0, step = 0,
lwd = 0.4, add = TRUE)

zoom out even more to give perspective on location:

juneau.out <- getNOAA.bathy(lon1 = -179, lon2 = -125,
lat1 = 75, lat2 = 48, resolution = 8)

plot(juneau.out, image = TRUE, land = TRUE, lwd = 0.1,
bpal = list(c(0, max(juneau.out), "grey"),
c(min(juneau.out),0,blues)))
# Making the coastline more visible
plot(juneau.out, deep = 0, shallow = 0, step = 0,
lwd = 0.4, add = TRUE)

That doesn’t look as good as the DEM, but hopefully it can give us the bathy transect.

# Example 1. get.transect(), without use of locator()
    get.transect(juneau, -134.8213, 58.4848, -134.7925, 58.4946) -> test ; plot(test[,3]~test[,2],type="l")
    get.transect(juneau, -134.8213, 58.4848, -134.7925, 58.4946, distance=TRUE) -> test ; plot(test[,4]~test[,3],type="l")

# Example 2. get.transect(), without use of locator(); pretty plot
    par(mfrow=c(2,1),mai=c(1.2, 1, 0.1, 0.1))
    plot(juneau, deep=-800, shallow=-10, step=1000, lwd=0.5, col="grey50",drawlabels=TRUE)
    #lines(nw.atlantic.coast)

    get.transect(juneau, -134.8213, 58.4848, -134.7925, 58.4946, loc=FALSE, dis=TRUE) -> test
    points(test$lon,test$lat,type="l",col="blue",lwd=2,lty=2)
    plotProfile(test)

# Example 3. get.transect(), with use of locator(); pretty plot
## Not run: 
    par(mfrow=c(2,1),mai=c(1.2, 1, 0.1, 0.1))
    plot(juneau, deep=-800, shallow=-10, step=1000, lwd=0.5, col="grey50",drawlabels=TRUE)
    lines(juneau)

    get.transect(juneau, loc=TRUE, dis=TRUE, col=2, lty=2) -> test
    plotProfile(test)

bathy.df <- highest.res %>%
  rename(longitude = x, latitude = y, depth = value)

bathy.df %>%
  write_csv("amalgaBathyDf.csv")

read.bathy("amalgaBathyDF.csv", header = T)

# Example 1. get.transect(), without use of locator()
    get.transect(highest.res, -134.8213, 58.4848, -134.7925, 58.4946) -> test ; plot(test[,3]~test[,2],type="l")
    get.transect(highest.res, -134.8213, 58.4848, -134.7925, 58.4946, distance=TRUE) -> test ; plot(test[,4]~test[,3],type="l")

# Example 2. get.transect(), without use of locator(); pretty plot
    par(mfrow=c(2,1),mai=c(1.2, 1, 0.1, 0.1))
    plot(test_spdf, deep=-800, shallow=-10, step=1000, lwd=0.5, col="grey50",drawlabels=TRUE)
    #lines(nw.atlantic.coast)

    get.transect(test_spdf, -134.8213, 58.4848, -134.7925, 58.4946, loc=FALSE, dis=TRUE) -> test
    points(test$lon,test$lat,type="l",col="blue",lwd=2,lty=2)
    plotProfile(test)

# Example 3. get.transect(), with use of locator(); pretty plot
## Not run: 
    par(mfrow=c(2,1),mai=c(1.2, 1, 0.1, 0.1))
    plot(test_spdf, deep=-800, shallow=-10, step=1000, lwd=0.5, col="grey50",drawlabels=TRUE)
    lines(test_spdf)

Also this: https://stackoverflow.com/questions/47047623/projectraster-raster-projection-of-bathymetry-data-noaa-nc-in-the-pacific